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I. INTRODUCTION 


Computational fluid dynamics has matured significantly 
within the past decade because of the development of 
increased computational capabilities and powerful 
computational techniques. Current problems being addressed 
include predicting flow separation over airfoils and 
post-stall flight characteristics. These areas are of 
interest because studies [Ref. 1] indicate increased lift 
and thus sustained flight are attainable when an airfoil is 
dynamically stalled, that is, its angle of attack is pitched 
to a post stall angle of attack rather than being initially 
placed at that high lift. 

Computational methods utilizing the full Navier-Stokes 
(N-S) equations are capable of addressing these issues, as 
are methods that include approximations to the Navier-Stokes 
equations. One method, the Interactive Boundary Layer (IBL) 
technique, developed by Tuncer Cebeci at Douglas Aircraft 
Company and at the California State University [Ref. 2], 
divides the flow over an airfoil into a viscous inner 
boundary layer and an inviscid outer layer. The 
characteristics of the inner flow are obtained from a 
numerical solution of Prandtl's boundary layer equation and 
the outer flow's characteristics are determined from Hess 


and Smith's panel method, and Fourier analysis and conformal 


mapping. The inner and outer layers are then redetermined 
by an interaction model that iterates between the two 
regions and marches downstream until the flow conditions 
have been satisfied at the boundary for both regions. The 
Cebeci IBL code uses Michel's criterion to predict 
transition from laminar to turbulent flow or transition may 
be prescribed. An algebraic (Cebeci-Smith) turbulence model 
is used. 

A full Navier-Stokes code developed by N.L. Sankar and 
his associates at the Georgia Institute of Technology [Ref. 
3] uses an implicit finite-difference procedure to solve the 
2-D Reynolds-averaged compressible Navier-Stokes equations 
in strong conservative form. The time-marching algorithm 
used is an Alternating Direction Implicit (ADI) procedure 
developed by Beam and Warming [Ref. 4] and implemented by 
Steger [Ref. 5]. The Sankar N-S code uses a body-fitted 
C-grid system and an algebraic (Baldwin-Lomax) turbulence 
model. 

The Cebeci IBL and the Sankar N-S codes are designed for 
different purposes. A low Reynolds number flow over an 
airfoil tends to be laminar until separation. The flow then 
transitions to turbulent flow and reattaches as turbulent 
flow. The Cebeci IBL code models this separation bubble if 
transition is specified within the separation bubble. 
Velocity profiles and skin coefficients are extremely 


important in analyzing these low Reynolds flows. Cebeci has 


developed codes for compressible oscillating airfoils, 
however, this Cebeci IBL code was developed for 
incompressible steady state flow only and thus does not 
predict the effects of unsteady flow nor compressibility. 

The Sankar N-S code was developed to address dynamic 
stall and its implication of increased lift. Therefore, the 
values of interest to date have been coefficients of 
pressure, lift, moment, and the effect of hysteresis on 
these values. However, the Sankar N-S code assumes the flow 
is fully turbulent, and therefore does not account for 
transition from laminar to turbulent flow. 

Neither dynamic stall nor transition within a separation 
bubble are easily quantified experimentally. Transition is 
a boundary layer phenomenon and the velocity profiles and 
skin frictions within the boundary layer must be measured to 
assure correct interpretation of the flow under 
investigation; surface pressures are not sufficient to 
accurately locate flow separation and reattachment [Ref. 6]. 
This is a time consuming, expensive process prone to error. 
Experimental methods include laser anemometry and hot wire 
probes [Ref. 7]. Disturbance of the boundary layer flow due 
to probes is undesirable and hot wire probes are normally 
unable to determine flow direction; therefore laser 
anemometry, although expensive and tedious, is increasingly 


being used. 


Dynamic stall is difficult to characterize due to the 
transitory nature of the phenomenon. Experimental 
techniques and apparatus include pressure transducers, hot 
wire probes, and laser doppler velocimetry [Ref. 8]. Flow 
visualization is also a very effective tool for studying 
both dynamic stall and separation bubbles [Refs. 9,10]. 

A good review of the current state-of-the-art 
computational and experimental aspects of aerodynamic flows 
is given in the proceedings of three symposia on this topic 


edited by T. Cebeci [Ref. 11]. 


It. OBJECTIVES 


The intent of this study is two-fold: to become 
familiar with computational fluid dynamic methods and to 
evaluate two codes to determine their range of 
applicability. 

Computational fluid dynamics consists of various 
mathematical methods and implementation schemes. A 
Significant portion of eralyede inherent in computational 
codes is empirical; therefore, the assumptions used strongly 
influence the results. It is important, when attempting to 
choose a computational code for a specific purpose, to be 
familiar with the significance of the analytical methods, 
assumptions made, and empirical models. Each code is 
different in these respects and must_ be analyzed 
individually and in detail to assure reliable, accurate 
results, especially when extending the flow regime or 
airfoil to conditions whose features are unknown. 

The two codes chosen, the Cebeci IBL code and the Sankar 
N-S code, are a good representation of two powerful, 
accurate methods that differ widely in computational 
approach. The Cebeci IBL code has been extensively tested 
in a variety of steady state conditions [Ref. 12] and the 
Sankar N-S code has compared well to experimental data for a 


pitching airfoil [Ref. 13]. However, the lack of steady 


state boundary layer data for the Sankar N-S code indicated 
that a more in-depth analysis of the applicability of the 
code was required. The Cebeci IBL code was used as the 
reference for the Sankar code. 
The analysis included the following: 
1. Assess C) for Reynolds numbers of 1.5 and 6 million 
2. Assess C, and Cr for a range of Reynolds numbers (1-15 
million at O degrees angle of attack) and angles of 
attack (0, 2, 4, and 6 degrees for 1.5M Re number and 
O, 4, 8, and 12 degrees for 6M Re number) 


3. Assess velocity profiles for Reynolds numbers of 1, 6, 
and 15 million 


4. Assess the influence of dissipation factors and 
grid size on the results of the Sankar N-S code. 


IIIT. MATHEMATICAL FORMULATION 


A. GOVERNING EQUATIONS 
Flow over an airfoil can be described by the velocity 
a fa + A Sa aie the pressure, the density, and the 
temperature. These six variables (u, v, WwW, p, 0, and T) are 
fully described by the continuity equation, the equation of 
state, p = oRT; the energy equation, 6Q - 6W = SE; and the 
three equations of motion. [Ref. 14] 
The continuity equation states mass is conserved; 
i.e., the flux of mass through a cube per time is equal to 
the time rate of change of mass. This is shown in Figure 


3.1 for the flow through the faces perpendicular to the x 


axis. Mathematically this is expressed as 


oew) Ax] AyAz - ae AY] AZAX - po sew) Az] AxAy 





(pAxAyAz) (32) 


Since the control volume is- fixed, AxAyAz is 


independent of time; therefore 


d(pu) , dev) , dfpw) _ 7: 


x + By + Oz. (3.2) 





or 


o u Ay Ax 


A (ou +2182) 0x] aydz 
y 





Net flux through face perpendicular to x-axis 


2 four Ax} AyAz 


Source: [Ref. 14:p. 106] 


Figure 3.1 The Flux of Mass Through a Cube 


92+ A+ pq@e0. (Gee) 


Alternately, 
Bet oe + ae + oe) = 0 54) 
ee 
ne + o(A-q) = 0 (3.5) 


The three equations of motion, one for each axis of the 
Cartesian coordinate system, are described by Newton's 
second law, AF = A(ma). The summation of the x components 


of surface forces on the element shown in Figure 3.2 is 





: a Du | GB) 
eae — eae "= (pAxAyAZ) pF 
with 
a0, 
AF =- 0, AyAz +: (cL, 1 aes Ax) AyAz 
OT 
- 1 Axdz + (T + I> Ay) bxhz 
oT Oy 
- t_. AxAy + (TL. + Az) AxAy . (3.7) 


STRESS-STRAIN RELATIONS 





Figure 3.2 X-components of Surface Forces on an Element 
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Dividing by the volume of the element yields 


Similarly, 


and 


mY, 
dy 


+ 


For Newtonian fluids 





for the y and z 











(3:8) 


foie) 


(31.0) 


with a single viscosity coefficient, 


the normal and tangential shear stresses are as follows 


[Ref. 15]: 


ig 
XY 


2u 


2 


Za 


zy 


(Sele) 


(echgullliey) 


(sel ley 


(Shape lls 


(eke) 


(S115) 


Substituting these stress definitions into the equations of 
motion, and assuming a constant viscosity corresponding to 
the mean temperature of the fluid ultimately yields the 


Navier-Stokes equations: 


2 2 2 
gp du du du LO re .4y fT Du 
Sax 7 Mt ee Goa ane aides 
OX dy dZ 
2 2. 2 a 
- <2 i i + = Pym +2 2 ty-q] = 0 = (3. 125) 
- 9x dy Zz Y 
Z 2 2 rn 
agp dW dW dW ua . _ Dw 
Se re i or ee oe 
aX dy dz 
or in vector format 
-Up + uv*g + Bv(v-q) = p <4 + p(q-V)q (3. 1) 
[Ref. 14] 


B. REYNOLDS STRESSES 

The Navier-Stokes equations are valid for laminar and 
turbulent flow. However, the complexity of turbulence has 
made it impossible to relate the motion of the fluid to the 
boundary | conditions and obtain an exact solution. 
Therefore, the turbulence must currently be modeled. OF 
Reynolds divided the turbulent flow into a mean motion and 


fluctuating, or eddying, motion as follows: 


12 


i) = el 
vy =a 
w=wWw 
pP=p 
p =p 
T=fT 


where the barred terms are the 


and the slashed terms are the 


the time averages of all 


fluctuations are equal to Zero: 


u! s., 14a) 


sa (3.14b) 
w' (2iliAie) 
io) (3.14d) 
oy (3.14e) 
dhs (3.14f) 


time-average of the component 
fluctuations. By definition, 


quantities describing the 


Rules for operating on mean time-averages are given 


below. F and g are dependent variables, and s is the 


independent variable x, y, Zz, or t. 
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(Ref. 16] 


The stresses caused by the fluctuations can be 
determined using the momentum theorem. consiaee an area dA 
with dA * pu * dat being the mass of incompressible fluid 
passing through the element in time dt. MThus, the flux of 


momentum in the x direction is 


dJ, = dA * pu? *dt; (3.15a) 
Correspondingly, 

aJy = dA * ouv * at (3:..15)7) 
and 

aj, = dA * opuw * dt. (3. Ses 


Calculating the time averages for the fluxes of momentum 


per unit time yields: 


a, = dAp u? (3.16a) 
dJy = dA p uv (3.16b) 
adJ, = dA op uw. (3.16c) 
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Utilizing the definition of turbulent flow and the previous 


rules yields 


dy, = ca - eae u'*) (3.17a) 
qa. =@A- p(u’v + u'v') (3.17b) 
ay, = dA - o(uew + u'w') (3.17c) 


Dividing these rates of change of momentum by area dA, we 
obtain stresses. The equal and opposite stresses exerted on 
the area by the surroundings are a normal stress, -(U@ + 
u'*), and two shearing stresses, -(uv + U'v") and -(Uw + 
u'w'). Thus, the superposition of fluctuations on the mean 


motion gives rise to three additional stresses 


Oo.=-pu', T' =-pu'v', Tt’ =-pu'w' . (3.18) 


The total stress tensor due to the turbulent velocity 


components of the flow are 


oO Tey i Vv u'w' 

Tey o — = -0 ie viw' ‘ (3.19) 
ci 13 of v'iw' ie 

SZ yz Z 





tS) 


The presence of fluctuations presents itself as an 
apparent increase in stresses (viscosity). These 
additional stresses over the mean or laminar stresses are 


termed apparent stresses or Reynolds stresses. 


B. TURBULENCE MODELING 

The presence of Reynolds stresses in turbulent flow 
introduces additional unknowns in the Navier-Stokes 
equations. Therefore, the Navier-Stokes equations, 
continuity, the perfect gas law, and the energy equation are 
no longer sufficient to completely define a solution. fThis 
is known as the closure problem and is usually resolved by 
turbulence modeling. [Ref. 17] 

A common method used is to relate the turbulent stress 
to the mean flow properties through empirically based 
algebraic formulas. An eddy viscosity, 4%, is defined in 
the same form as the laminar viscosity. Previous models 
related surface boundary conditions to points in the fluid 
away from the boundaries through wall functions. This 
avoided modeling the direct influence of the eddy viscosity; 
however, it is only applicable in regions where the Reynolds 
number is high enough for viscous effects to be unimportant 
or where universal wall functions are well established. [In 
turbulent boundary layers at ise Reynolds numbers, in 
unsteady or in separated flows, or in three-dimensional 
flows, the flow close to the wall must be described. [Ref. 


18] 
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A common algebraic turbulence model divides the flow 
into an inner and eunee layer. The inner layer is defined 
by a modified mixing length formula that utilizes some 
damping function. The outer layer includes the wake and 
another damping function. [(Ref. 17] 

Other empirical methods currently in use, such as the e? 
method, predict transition. The el method is a stability 
method, based on linear stability theory. It assumes 
transition begins when a small disturbance is introduced at 
or below a critical Reynolds number. The transition is 
amplified by e”%. This method allows greater generality of 
the flow, however the formulation still relies on empirical 
terms. [Ref. 19] 

The accuracy of turbulence models are limited by the 
accuracy of the empirical constants. Caution must be taken 
when using a model under different conditions, i.e., a 
different flight regime or a radically different airfoil. 
The turbulence models mentioned above can be fairly simple; 
a more complex model is still not generally usable because 
of the computation costs involved and the uncertainty of the 
constants. 

The influence of turbulence and the transition from 
laminar to turbulent flow on the airfoil need to be 
understood and accurately modeled for a good description of 
the flow over the airfoil to be detailed. As experimental 


methods continue to improve and as computational methods 
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utilize the data, improvement in the detail of the flow 
field and in flow prediction will follow. 
1. Cebeci-Smith Turbulence Models 

The turbulence model used in the Cebeci Interactive 
Boundary Layer (IBL) is a simple algebraic eddy viscosity 
expression. Simple algebraic models seem to adequately 
predict turbulent flow for wall boundary layer flows in 
which the Reynolds shear stress and frequency do not change 
rapidly. However, if the rate of change of shear stress or 
the frequency is large, turbulence models are not currently 
satisfactory. [Ref. 2] 

The Cebeci IBL code utilizes the algebraic 
eddy-viscosity formulation of Cebeci and Smith [Ref. 20]. 
The turbulent eddy viscosity, Vy, for wall boundary flows is 
defined by two separate formulas; one for the inner region, 
based on the Van Driest approach, and the other for the 


outer region, based on a velocity defect approach: 


2 ,du 
{0.4y[1 - exp(-y/A)]}" |a7ly,, for Os ¥<syY, 


v, = : (3.20) 
i Cee ee. 7 for gy YS 3 
0 
where: 
Ly 
du 
A 26v| (v =) and 
ee: re 
= 5: 
bao) 
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The continuity of the eddy viscosity defines y,; the 
expression for the inner region is used outward from the 
wall until it agrees with the outer region, which is then 
used. Ytyr 1S an intermittency factor which allows for a 
transition region when progressing from laminar to turbulent 


flow. It is given by 


w -1.34 x 
Y = 1 - exp[- —=-.R (x - xX ) f de, (3221) 
jeg 2 teTs u ; 
G Vv tals x e 
Yer tr 
where the transitional Re number, Re = (Ugx/Vty)- The 


empirical constant Gy+; is dependent on the Reynolds number 
of the flow. High Reynolds numbers flows indicate Gyty = 
1200, lower Reynolds number flows seem to be better modeled 
by lower values of Gey;y-. [Ref. 12] 

The parameter q in the outer region is given by 
®© = 0.0163/F2°> where F is the ratio of the normal stress 
turbulent energy to the shear stress turbulent energy 
evaluated at the point of maximum shear stress. This can be 


expressed as 


pees 


= = )(ur? - vt)au/ax | (26) 


- u'v' du/day eee ae 


The ratio of the time-averaged quantities are 
assumed to be a function of Rp = 1,/(-u'v') max which can be 


represented by 
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6 
———— <a: 0 
( oe - oo 1+2R, (2-R,) Rp 
aa | 2 Me | = (8523) 
- u'v' —— 2Rn 
(on) 
max TR, Rn > 1.0 
Thus, the expression for a is 
i 0.0168 ; [Ref. 2] 
[1 - 8 (3u/ax)/(du/ay)]°"> 
(3324) 


Note that the value of Y has the effect of reducing 
the eddy viscosity away from the airfoil surface. This 
turbulence model does not take into account the wake region, 
nor is it validated for separated flow. 

2. Baldwin-Lomax Turbulence Model 

The Sankar Navier-Stokes (N-S) code also uses an 
algebraic eddy viscosity model, the Balwin-Lomax Turbulence 
Model [Ref. 21]. It is based on the Cebeci-Smith two layer 
model [Ref. 22] used in the Cebeci IBL code and may be 


expressed as 


2 
{.4y[1 - exp(-y/A)]} |u| for O<y<y. 
= eZ 
Ve (SEZ) 
BO RUSE LG's) aes Feb LY? for yi<y <6 
where 


A = 26. 
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The inner region is the same mixing length formula 
of the Cebeci-Smith model, simplified. No intermittancy 
factor is included (flow is calculated as wholly turbulent), 
A is a constant rather than being dependent on viscosity and 
velocity gradients, and the velocity profile (du/dy) is 
replaced with the product of vorticity and density. 
(Ref. 13] 

The outer region is based on the wake function, 


Fwake and the Klebanoff intermittancy factor, Frjep(y).- 


PRet. 22) 
: 2 
Bee min (y Pigee 23 -V. Uaie/F a (S220) 
6 -l 
Feb ¥) = (1+ 5.5 (.3y/y J (=27) 


The quantities yay and Fmay are the maximum values obtained 


from the function 


F(y) = y |w| {1 - exp (-y'/a’)} (3.28) 


which is a form of the mixing length formula used in the 
inner region. 

Ugif is the difference between the maximum and 
minimum velocity of the velocity profile. Fyarpeis similar 
to the y of the Cebeci-Smith turbulence model and thus also 


reduces eddy viscosity away from the airfoil surface. 


Pa 


This model has been used in separated flows and in 


the wake, however its validity is not assured in these regions. 


ODP 


IV. THE BOUNDARY LAYER EQUATION 


Prandtl clarified the influence of viscosity in high 
Reynolds flows by simplifying the Navier-Stokes equations to 
yield approximate solutions. He divided the air flow over a 
body into two regions: 


1) The region near the surface where viscous forces 
dominate. 


2) The rest of the flow where inertia forces dominate; 
this region may be considered frictionless and 
potential. 

Consider a 2-D incompressible flow over a body. Most of 
the flow is moving at free stream velocity. However, at the 
surface the velocity is zero, increasing to free stream at 
some distance from the surface as shown in Figure 4.1. [In 
this first region, called the boundary layer, the velocity 


gradient normal to the wall, o/sy, is very large, as is the 


shearing stress, 
t=yp 2 (4.1) 


The two regions are not distinct, but are usually divided at 
the streamline where the velocity reaches 99% of the free 
stream velocity. 

Simplification of the Navier-Stokes equations will be 
accomplished by doing an order of magnitude analysis of each 


term. The following assumptions are made: 
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Pressure decrease Pressure increase 


Figure 4.1 Velocity Profile 
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1) A flat wall coinciding with the x-direction, with the 
y axis perpendicular. 


2) Ones 
Sj be Al 


where 6 is boundary layer thickness and L is a linear 
dimension of the body so selected to ensure du/ox = 1 
under the region in consideration. 


3) The Reynolds number, R = ULp/ w= is large. [Ref. 
16] 


The Navier-Stokes equations are rewritten in 
dimensionless form: 


- Velocities are with respect to the free stream velocity, 
u. 


- Linear dimensions are with respect to a characteristic 
length, L. 


- Pressure is divided by pU?. 
Retaining the same symbols for the dimensionless quantities 
as for their dimensional counterparts, and writing the order 


of magnitude under each term yields for continuity: 


Chel Tr 2A 4.2 
ny Oy 0 (4.2) 
Jk 1 


and for the Navier-Stokes equations: 


Z eZ 
du du co) ae Cee ae | 
Gir x uxmtv—aHr at alst+ —) (4.3a) 
ox dy ox Rony? a ; 
Ll 6% a, 
6 
dV dV oe! 3¢v 3¢y b 
; oe ee ee 4.3 
cule yrs 1s 5 a sp Gell gp oP oy) ( ) 
-) oy 
ions f $ 6° 8 = 


At the outer edge of the boundary layer u equals the 


steady flow U(x). The viscous terms no longer dominate and 


thus, for the outer flow 


u@--52 (dimensional) (4.4) 
or rewritten in the form of Bernoulli's equation 
D sey Zane ur = constant - (4.5) 
Returning again to dimensional quantities, the 


simplified Navier-Stokes equations, known as Prandtl's 


Boundary-Layer Equations may be written: 


ou. OV _ 

a os © (4.6a) 
du ou 1 dp 3~u 

Max” Vay” pa V2 (4-6 


with the boundary conditions y = 0: u= 0 v= 0; y = inf: 


u = U(x). Also, a velocity profile at the initial section, 


xX = Xo, must be prescribed. [Ref. 16] 
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V. THE SANKAR NAVIER-STOKES METHOD 


The Navier-Stokes Code utilizes the governing equations 
in conservative form, a body fitted coordinate system, and 
an Alternating Direction Implicit (ADI) procedure [Ref. 23]. 
The governing equations in conservative form have the 
coefficients of the derivative terms constant. The 
conservative form ensures the conservation of the physical 


properties. 


A. GRID GENERATION 

The requirements of a grid in the physical plane and in 
the computational plane are conflicting--therefore a grid 
transformation is advantageous. For ease in computation, 
equal spacing of grid points is desirable; however, the 
physical grid needs to be clustered so that the boundary 
layer and sharply curving surfaces such as the leading edge 
contain enough points so as to be adequately defined. The 
boundary conditions must be accurate and should be contained 
on rectangular surfaces in the computational plane. Also, 
the grid should be smooth with few discontinuities. 

The present code uses an algebraic C-grid which 
generates a sheared parabolic coordinate system [Ref. 23], 


first proposed by Jameson [Ref. 24]. 
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First, two points, T and N, are defined on the desired 
airfoil in the X-Z plane as shown in Figure 5.1. These 


points are defined by the complex values 


Zp = Xp + lyr (5.1la) 


2N xy + lyn (5. ib 
respectively. T is located at the trailing edge of the 
airfoil and N is located half way between the leading edge 
and the center of curvature of the leading edge. 

Next, the airfoil in the physical plane is transformed 
as shown in Figure 5.2. A trailing edge vortex sheet shape 
is assumed to leave the trailing edge smoothly by running 
tangent to the mean camber line at the trailing edge. The 
airfoil and wake are then mapped onto the plane by using the 


following transformation, 


G =vZ-~ Z2y .- (5.2) 


The NACA 0012 airfoil shape transforms to that shown in 
Figure 5.3. Cubic interpolation defines additional points 
to smooth the surface. 

The far field boundaries are mapped in the physical and ~ 


computational planes as shown in Figure 5.4. 
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Figure 5.1 Defined Points on the Airfoil Surface 
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Source: (Rete isi pease 


Figure 5.2 Symmetric Airfoil in the Physical Plane 
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Source: ORefs I3:p. 19) 


Figure 5.3 Symmetric Airfoil Unwrapped in 
the Intermediate Plane 
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Figure 5.4 Far Field Boundaries 
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A sheared Cartesian coordinate system is then 
constructed in the zz plane. It consists of straight 
vertical lines which then contain specified clustered 
spacings defined by a stretching function. This allows the 
grid size to increase normal to the surface. 

Lastly, the grid is mapped back to the physical plane. 


The points on the computational plane, ¢, are given by 
Cc = €+ in (5.3) 
and on the physical plane by 
X- xX, = & + 7 (5.4a) 


Y= Yn = 2a (5.4b) 
(Refs. 13,23]. 

The present method uses a grid containing 161 points in 
the € direction and 41 points in the Ndirection. The final 
grid in the physical plane is shown in Figure 5.5 and a 
detail of the grid around a NACA 0012 airfoil is presented 


in Figure 5.6. 


B. INITIAL CONDITIONS AND BOUNDARY CONDITIONS 

The initial conditions for viscous flows are the free 
stream conditions. Viscous dissipation inherent in the 
equations will minimize seeers in this approximation after a 
sufficient time. Inviscid flows require the proper 
combination of "artificial" dissipation and boundary 


conditions for the correct result. 
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Figure 5.5 Final Grid in the Physical Plane 
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Boundary conditions must be defined on the airfoil, at 
the far field, and in the wake. On the solid surface, two 
conditions must be satisfied; no penetration of the surface 
and the solid and fluid have the same velocity at their 
boundary (the no slip condition). Adiabatic flow is 
assumed. 

In the far field, the flow is represented by the linear 


small disturbance equation 
(1 - M,)?7%&y + %yy = 0 (eS) 


where ? is the perturbation potential and x and y are the 
physical plane coordinates. This model is used instead of 
specifying flow conditions at infinity to compensate for the 
loss of lift experienced when the boundary is not placed far 
enough away from the solid surface [Ref. 23]. 

In the wake, the grid procedure produces a "cut" along 
the coordinate line from the trailing edge to the far field 
boundary. Here, the flow properties are averaged from above 


and below. [(Refs. 3,13] 


C. NUMERICAL FORMULATION 

The coupled and non-linear governing equations are 
solved using an alternating direction implicit time marching 
procedure similar to that developed by Beam and Warming 
[Ref. 4] and as used by Steger [Ref. 5]. The governing 


equations are assumed to have a solution at some time t, and 
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the solution is then advanced to some time t,4 1 using Euler 
implicit time differencing. The equation is then linearized 
using Taylor expansion and reduced to ae series of 
one-dimensional problems using the Beam and Warming 
approximate factorization method [Ref. 5]. 

The procedure is not fully implicit since the viscous 
terms are lagged by one time step. Artificial dissipation 
is included in both second and fourth order terms to obtain 
accurate surface pressures. The ADI procedure is limited in 
time step because of accuracy but the explicit boundary 
conditions further limit the step size due to stability 
considerations. 

1. Governing Equations 

The governing equation in the computational plane is 


given as 


pq +o, E+ oF = Re“1(9,R + an S) (5.6) 
where 
q=q/J (5.7) 
E= (4G + &f + EyF)/J (5.8) 
F = (ma + nyE + nyF)/J (5.9) 
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Tyx 


Txy = 


eyey 


a 


a (E,R + EyS/J 


4) 
It 


(nyR + nyS)/T 

Cee) ( gx, + Hy, ) zi \ (EyV pt Nyv >) 
uC (éyU,tnyu, ) + C&vet hv) ] 

(At+2 y) (Eyv + ryv) Te M(E xu, + nyu, ) 

+ Viyy + kPr71(y-1) (8,3, a2 y 3, a2) 


+ Viyy + kPr71(y-1) (by 95 a2in y 9,a2) 


J is the transformation Jacobian: 


T= Exenym Gy = (%Y_ Ye) 


Also, 


= J 
Ex Yy 
Nx = ne es 


SS) 


(5.10) 


(5.11) 


(5.12a) 


(5.12b) 


(5.12c) 


(5.13a) 


(5.13b) 


(5.14) 


(5.15a) 


(5.15b) 


(5.15c) 


‘an 


= “x Ex y by (5.15e) 


Details of the derivation of the governing equations in the 
physical plane and transformation of the governing equations 
to the computational plane are contained in References 13 
and 23. 
2. Time Differencing and Linearization 

The governing equations cannot be solved directly in 
the form of Equation (5.6) because of their non-linearity. 
This is overcome by writing the flow quantities op, pu, Pv, 
and e at their new time level as their value at the known 


time level and their increment, i.e., 


n+l 1 


5 So Noe (5.16) 


The variable q can therefore be written at the new 


time level tp41 as 


or 


So 


Aqnt1 = qntl a qh ’ (5.17b) 


A Taylor series expansion of q! backward about qntl yields 
gM = gMtl —- atseqgnt? + 0( At?) (5.18a) 
qntl -q = At dqntl i O( At?) a Agnt1 (5.18b) 


Using the same procedure with central differencing for the 
Spatial derivatives and substituting into the above equation 


yields 


Aqntli = At (6, Entt + Seer) nm dRe~1(5,R0+1 


+ as gntl) + 0 (At) (5.19) 


where Or and 6, are second order accurate difference 
operators. 

Backward differencing is first order accurate and 
central differencing is second order accurate, therefore the 
Equation (5.19) is first order accurate in time and second 
order accurate in space. 

The governing equation, Equation (5.6), is still 
non-linear. So that it can be solved directly rather than 
iteratively, it is linearized. First E and F are rewritten 


using a local Taylor expansion at q?: 
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ENt+l = EN + [oqE)MaqMt! + 0( At?) (5.20a) 
FN+L = PN + [aqF]MagMtt + o( at?) (5.20b) 


where [dqE] and [9qF] are Jacobian matrices. (Ref. 23] 
Substituting Equation (5.20) into the time differenced 


governing equation, Equation (5.19), yields 


({I] + Até,[ ae)" + at 6 (aaF]™) ( aayn*? 


= -At(5, E45 FM) + atRe™1(6, oe HSU) (5.21) 


E 


This can be expressed as 


({T]+at 6A] + at 6 (B]) (Aayht> = (R}M (282) 
where: 

[A] = (3 qE)” (Se22) 

([B]) = (9 qF]” (5.24) 

{R}? = Sis eee) + AtRe~? (6 RATI+ 3 snr) (5.25) 


The governing equations for the unknown vector 


{Aq)"*1 are now linear and may be solved numerically. 
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3. The Alternating Direction Implicit Procedure 


An approximation factorization method, developed by 
Beam and Warming [Ref. 4] is now used to solve the governing 
equation Equation (5.22), which has been linearized in 
[Aq]mtl, Although Equation (5.22) could be solved directly, 
Beam and Warming's method reduces the large (and thus 
costly) matrix into a sequence of one-dimensional problems. 
The governing equation is approximated as 


({I] + Até_[A]) ({I] + At 6 (B]) (aq}nr* = {R}n (5.26) 


3 


which is then rewritten as two equations 
(bee I) a At 5,[A]) (Ag*} = {R}" (5.27a) 
({I] + At 6 (B]) (aaynr? = {aq*) (5.27b) 
Since Se and 6, are central difference operators, Equations 
(5.27) are systems of block tridiagonal matrix equations 
composed of 4x4 submatrices. These can be solved by one of 
several methods, such as LU decomposition. Once {Aq}"*1 is 


obtained, {q}"t1 can also be obtained from 


(qyMtl = (q)Ph + (aqpntt (5.28) 


42 


Boundary conditions at {Aq} as well as the defined vector 
{Aq*} must be defined. 
4. Stability and Accuracy Considerations 

The ADI approach is implicit with explicit boundary 
conditions and viscous terms that are lagged by one time 
step. Implicit numerical methods theoretically are 
unconditionally stable with the size of the time step 
limited by accuracy rather than stability. The time step 
stability limit imposed by the explicit boundary conditions 
must therefore be less than the accuracy time step limit. 
fRef. 4] 

A linear stability analysis for an explicit 
procedure is performed [(Ref. 23] to determine the maximum 
time step possible, then a very conservative estimate is 
used. When calculating steady-state flow, the convergence 
of the calculations can be improved by introducing a 
variable time step. This allows small time steps where the 
grid must be highly clustered, such as the boundary layer, 
in regions of shocks, and near stagnation points, and large 
time steps elsewhere where the grid is larger. [Ref. 23] 

Viscosity slows a flow down by dissipating energy. 
Mathematically, this results in a reduction of flow field 
gradients. The mathematical formulations used to include 
viscosity in a flow can therefore also be applied to 
Suppress errors inherently generated in certain methods. 


[Ref. 25] 
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The capability of including” this “artiuederal 
viscosity" as used by Murman and Cole [Ref. 26] is 
implemented through backward differencing in a Taylor series 
expansion. The truncation term mimics viscosity and is 
knewn as artificial viscosity. When the method includes 
this dissipative term it is known as implicit artificial 
viscosity. Often however, more dissipation is required for 
convergence or stability, or because it is advantageous to 
apply selective dissipation. Then, explicit artificial 
viscosity or explicit dissipation is included in the 
numerical formulation. 

Central differencing in a Taylor series expansion 
decouples the even from the odd terms, causing high 
frequency errors. The spatial derivatives are formulated 
through central differencing and thus contain these errors 
which influence the solution accuracy when using large time 
steps. The time derivatives are formulated through backward 
differencing; their artificial viscosity terms and the 
viscous terms suppress this error somewhat. To further 
dissipate the high frequency errors, especially in high 
Reynolds number and inviscid flows, both fourth order and 
second order artificial dissipation is explicitly applied to 
the right hand side of the descretized governing equations 
as was done by Jameson [Ref. 24]. [Ref. 23] 

The fourth order terms alone lead to overshoots in 


the vicinity of shocks. The second order dissipative terms 
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correct this but tend to smear results in regions such as 
the leading edge. This is resolved by implementing second 
order dissipation only in regions of high pressure 
gradients. Including artificial dissipation in the right 
hand side of Equation (5.26) also corrects for incorrect 
initial conditions after a sufficient time step. 

To allow the viscous terms to be explicitly 
modulated and to remove any explicit stability limits, 
artificial dissipation is also implicitly included in the 
left hand side of Equation (5.26). The inclusion of 
explicit and implicit artificial dissipation yields (Ref. 
23] 


Geet ACS. [A}- Eat 6-2) ({I]+ At 6 (B]- Egat Sane) { Aq}yntt 


E 


(Jq) (5.29) 


= ota a a oars fia 


where is a function of the maximum pressure gradient and D 
and D;_1 are either second or fourth order dissipation. The 
artificial dissipation terms do not affect the accuracy of 
the formulation since all terms are of the same order as or 
smaller than the truncation errors associated with the 


spatial and time difference formulas. 
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5. Application of Boundary Conditions 


All of the boundary conditions are explicitly 
applied after the ADI sweeps at each time step. The 
explicit method was used because of its ease of 
implementation in the code even though implicit boundary 
conditions are more desirable because of accuracy and 
stability considerations. [Ref. 23] 

On the airfoil surface the no-slip condition and the 


assumption of no flow through the surface correspond to 


Ip/an = 0 do/3n = 0 (5.30) 


Adiabatic flow is assumed. A two point extrapolation of the 


above yield p and op to be [Ref. 14] 


Pi, = (4EeS as Pi 3)/3 (5.31a) 


Pi, dupe AuDi | 2 eee oie (5.31b) 


Internal energy and the coefficient of pressure may be 
obtained from the calculation of p and ep. The incremental 
quantities {,Aq*} and {A qyntl are also assumed zero on the 
solid surface and are solved in a similar manner as for p 
and». The far field boundary conditions are assumed to be 


free stream plus the disturbances caused by the far field 
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not being placed at infinity as described in Section V.B. 


Internal energy is evaluated from 
2 - eee) (5.32) 


The downstream boundary conditions are extrapolated from the 
adjacent interior points. At the boundary, pressure is 
taken to be freestream. 

The cut inherent in the grid system divides the flow 
above and below the line emanating from the trailing edge. 
The flow quantities on the cut are obtained by averaging the 
values of the interior points above and below the cut. This 
is acceptable because of the denseness of the grid in this 


region. [Ref. 23] 


D. USE OF SANKAR!S N-S CODE 

Reference 13 details the use of Sankar's N-S code for 
both steady and unsteady flows. The code is submitted to 
the NASA X-MP Cray via Job Control Cards. JCL options are 
selected so to access files where output data is stored or 
to send data to a specific directory. Job Cards also 
provide account and time limit information. 

The main program contains all of the "Write" statements, 
which may or may not be required. Output information 
includes input data, airfoil coordinates in the physical and 
computational planes, grid information, residuals, pressure 


and skin friction profiles, coefficients of lift, drag and 
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moment, and velocity profile information. The time step 
intervals at which these are printed is specified within the 
program. 

The majority of the program inputs are located in data 
cards. These inputs and their definitions are decribed in 
Table 5.1. See Reference 13 for a detailed program 
description. The values used when implementing the Sankar 


N-S code are given in Table 5.2. 
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TABLE 5.1 


INPUT PARAMETERS FOR THE SANKAR N-S CODE 


INPUT 


PITCH 


PLUNGE 
FNU 
FNL 
ZSYN 
ae 


ae 


DESCRIPTION 
Number of x coordinate locations 
Number of y coordinate locations 
Size of time step 
Explicit artificial viscosity term 
Mean angle of attack 
Amplitude of oscillation 


Angle below which modified turbulence 
model is used 


Reduced frequency 

Mach number 

Number of time steps 

Reynolds number in millions 

Distance of first point off of the wall 
Time calculation flag 

Plotting file flag 


If set, stored values are read to 
continue iteration 


Set for dynamic stall, indicates change 
in AOA 


Set for up and down motion of airfoil 
Number of upper airfoil coordinates 
Number of lower airfoil coordinates 
Flag for symmetric airfoil 

X airfoil coordinates 


Y airfoil coordinates 
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TABLE 5.2 


VALUES USED FOR SANKAR STEADY N-S CODE 


INPUT VALUE 
IMAX L57 

KMAX 41 

DT 1.00 

WwW 10-12 
ALFA 0-12 
ALFA1 0 

ALFAIL ¢) 

REDFRE 0.0 

AMINF ri2). Soe 
FNSTP 2000-8000 
REYREF => 
DNMIN -000005-.0001 
TSTART 0 

RSTRT 0 

PITCH FALSE 
PLUNGE FALSE 

FNU 33 

FNL 33 

ZSYN 0 

XP X airfoil coordinates for NACA 


0012 airfoil 


YP Y airfoil coordinates for NACA 
0012 airfoil 
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VI. THE CEBECI INTERACTIVE BOUNDARY LAYER METHOD 


The Cebeci Interactive Boundary Layer (IBL) code uses a 
two layer approach; an outer inviscid layer and a viscous 
inner boundary layer. In this it follows Prandtl's boundary 
layer theory by assuming the inner viscous flow is the only 
region where viscous effects are important. The remaining 
outer region is dominated by inertia terms and can be 
assumed inviscid. The inviscid outer flow characteristics 
are determined through Hess and Smith's panel method. The 
flow is assumed to have a vortex and source distribution 
such that it gives correct circulation and velocity over the 
airfoil. The airfoil surface is replaced by a distribution 
of panels that satisfies the Kutta condition. 

The inner viscous flow utilizes the boundary layer 
method to determine the flow characteristics. A direct 
boundary layer method is used in regions of large viscous 
stresses such as near the leading edge. This method 
prescribes an external velocity and requires the no slip 
condition to be satisfied. In the rest of the flow, an 
interactive boundary layer method is used. Here, the edge 
boundary conditions prescribe a combination of displacement 
thickness and external velocity. An iterative technique is 


used to solve for this flow. [Ref. 27] 


ag 


The inner and outer flow are coupled through an 
interaction model. The boundary location and velocity are 
unknown and are solved simultaneously through an iteration 
between the inner and outer flow equations. 

Cebeci's IBL code calculates where transition from 
laminar to turbulent flow occurs and incorporates a 
smoothing function. It also predicts separation and allows 
laminar separation and then transition to turbulent flow and 


subsequent reattachment of the turbulent flow. 


A. VISCOUS INNER FLOW 
1. Direct Boundary Layer Method 

The direct boundary layer method can be used in 
regions of the boundary layer where the viscous effects have 
not yet strongly influenced the flow. This usually implies 
the stagnation point and the airfoil leading edge. The 
advantage of minimal viscosity influence is the capability 
of defining a stream function, ~, that satisfies the 


continuity equation 
Oe — cy oy and v= - d/ dx (6.1) 


The momentum equation is subjected to the 


Falknar-Skan transformation [{Ref. 2] 


€ (6.2a) 
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i 
VU, VX 


(x,y) (6.2b) 





£(x,n) = 


The normal coordinate y and the stream function w are scaled 
with respect to the external velocity for convenience and 
accuracy. The boundary conditions of no flow penetration 
through the wall and no slip condition at the wall are also 
transformed. The resulting momentum equation and boundary 


conditions are given as 


(bE"')" + S(mel) ££" + m{1-(£")7) = x(é" SE - gt) (6.3) 
where: 
be- Jot Ye/ V 
m= (x/u, ) (du, /dx) 
and 
n = 0: iC 0) = O:, £(x,0) = 0 


heater: £ (xX, ey — L 


Primes denote differentiation with respect to n This is a 
third order partial differential equation and cannot be 
solved directly. Therefore, the box method developed by 


Keller [Ref. 28] is used. [Refs. 12,30] 
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The box method reduces the governing equations to a 
first order system through introduction of two dependent 
variables. The flow properties are then evaluated only 
discretely by defining the solution domain as a rectangular 
mesh. Instead of solving continuous functions, all 
parameters are approximated in terms of nodal values and 
their location on the mesh. The domain of dependence is 
substantially reduced and the overall solution scheme 
simplified by solving for the nodal values through central 
differencing. 

The resulting nonlinear system is solved by Newton's 
method. This iterative procedure linearizes the variables 
at location i by rewriting the value at i as a sum of the 
value at location i-1 plus some incremental value. 
Substituting into the governing system of equations results 
in a linear system of the unknown incremental values which 
are repeatedly solved until they are small enough to be 
neglected. 

2. Interactive Boundary Layer Method 

Most of the airfoil is influenced by viscosity, thus 
the direct boundary layer method can no longer be used. 
Instead, an interactive boundary layer method is used. It 
is effective even in regions of rapid flow acceleration, 
boundary layer separation, and zero skin friction. Both the 
boundary layer displacement thickness and the external 


velocity are now unknown. The Mechul function method is 
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used to solve the flow under these conditions by writing the 
edge boundary condition as a sum of the inviscid velocity 
distribution and the perturbation velocity due to viscous 
effects. The perturbation velocity, Su,(x), is determined 


from the Hilbert integral (Refs. 12,27] 


“b (al do 
Gug(x) = >) ay (,6*) = (6.4) 
x 


a 


where: 


dige5”) . 

ao 18 the blowing velocity. 
The solution method follows the direct boundary layer method 
with several exceptions. The Falknar-Skan transformation, 
with its constant boundary layer thickness, is no longer 
applicable; nor is using Ue as a reference velocity since u, 


is now unknown. Instead, 





LS 
n = oe Y (6.5a) 
£(x,n) = —~ (x,y) (6.5b) 
ees 


where the reference velocity is now taken as u,, the free 
stream velocity. The solution is again a downstream 
marching technique. The FLARE (Flugge-Lotz and Reyhner) 
approximation is used to continue the integration through 


regions of backflow. In these regions, where the velocity 
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in relation to the forward velocities is assumed small, the 
streamwise convection term udu/ ox is set equal to zero. 
3. Transition Model 

For the IBL technique to be successful, the 
displacement thickness must be accurately determined. This 
is dependent on an accurate solution of the laminar and 
turbulent flow equations, the transition region between 
them, and when applicable, separated flows. 

The values of flow parameters associated with 
laminar and turbulent flow differ greatly: the boundary 
layer thickness, momentum thickness, skin friction, velocity 
profile, and drag are all influenced by § increased 
turbulence. For a code to accurately model the boundary 
layer flow, both regimes must be modeled as well as the 
transition between the two and separated flow. The 
influence of the ePaHaeIOH location and length of 
transition on the accuracy of the solution, especially for 
low Reynolds flows, has been demonstrated by Reference 12. 

The Cebeci IBL code uses Michel's criterion to predict 
onset of transition. Michel begins transition when the 
local Reynolds number based on the momentum thickness, Ros 
is related to the length Reynolds number by the empirical 


equation 





225 
Re = 1.174(1 + R. ) R (6.6) 
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where: 


Rp = u,6|v ¢ and 


ye) 
6 


The transition model is also highly empirical and is 
as previously given in Equation (3.20). 

As implied by the above discussion, turbulence onset 
and its generation mechanisms are not thoroughly understood. 
Nonetheless, it is evident that numerical methods need to 
include transition capabilities to begin to accurately model 


the boundary layer flow. 


B. INVISCID OUTER FLOW 

Hess and Smith developed a panel method where the flow 
is represented by a series of vortices and sources. (Ref. 
30] They assume the vortex strength to be constant and 
distributed over the surface such that the correct 
circulation results. The velocity field is then modeled 
through a source distribution that forces the velocity to be 
everywhere tangent to the surface (the no penetration 
condition). This method is simplified by defining nodes on 
the airfoil surface and connecting them with straight line 
panels. Obviously, greater accuracy is achieved by 
increasing the number of panels on the airfoil. [Ref. 12] 

The Kutta condition must also be satisfied. The Kutta 
condition insures a unique solution by imposing zero loading 


in the region of the trailing edge. The three basic 
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principles of the Kutta condition are that the circulation 
about the airfoil is such that the flow leaves the trailing 
edge smoothly, the trailing edge is a stagnation point if 
the trailing edge is finite, and the upper and lower 
trailing edge velocities are finite and equal if the 
trailing edge is cusped. The vortex strength is determined 
from the Kutta condition and the source strength can then be 


calculated from the vortex strength. [Ref. 30] 


C. STRONG INTERACTION MODEL 

The inner and outer flow influence each other and thus 
cannot be solved separately if viscosity effects on pressure 
are large. The strong interaction model couples the 
boundary layer and the external viscous flow by allowing 
both the displacement thickness and the pressure (which is a 
function of external velocity) to be unknown. An iterative 
Simultaneous solution is then achieved by alternating 
between the viscous and inviscid flow equations until 
convergence is achieved. 

The solution method is based on conformal mapping and 
Fourier analysis techniques [Ref. 31]. The airfoil is 
mapped onto a circle through a series of conformal mappings 
and application of the fast Fourier-transform algorithm. It 
then is mapped onto another circle that includes’ the 
boundary layer and so models the inviscid portion of the 
flow as that over an airfoil whose boundaries have been 


displaced by the viscous boundary layer. At the surface of 
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the circle, the no penetration condition does not apply. To 
account for this, a nonzero normal velocity (blowing 
velocity) is prescribed. [Ref. 29] 

The boundary layer is solved from the surface to the 
outer boundary. The outer boundary conditions are defined 


by the interaction law 
e e K 
Gaye = (uly eat ge ey Ci-l [6 #)K/K-l oy, Sx) kK, K-1 (6.7) 


The solution of this is an approximation that requires 
several sweeps over the upper and lower surfaces to achieve 
convergence. 

Convergence between the boundary of the two methods 
is checked and the procedure is repeated by updating the 
product of the external velocity and displacement thickness 


until desired convergence is achieved. 


D. USE OF CEBECI'S IBL CODE 

Reference 29 details the use of the Cebeci IBL code. 
The code is submitted to the Naval Postgraduate School IBM 
mainframe via Job Control Cards where account information, 
running time, and size commands are set. 

Output information includes input data, coefficients of 
lift and drag, airfoil coordinates, shear stress, skin 
friction, displacement and momentum thickness, and velocity 
profile information. The inputs to the program are located 


in data cards and are defined in Table 6.1. See Reference 
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29 for a detailed program description. An example of input 


parameters used is given in Table 6.2. 


60 


TABLE 6.1 


INPUT PARAMETERS FOR THE CEBECI IBL CODE 


INPUT 


IPT (1) 


IPT (2) 


IRSTRT 
IGLMAX 
IPRINT 
RL 
XTRL 
XTRU 
ALPO 


GTR 


MPTS 
XP 


YP 


DESCRIPTION 


Flag for specification of lower surface 


transition 


Flag for specification of upper surface 


transition 

Flag for restarting solution 
Number of sweeps 

Flag for output printed 

Reynolds number in millions 

Lower surface specified transition 
Upper surface specified transition 
Angle of attack 


Empirical constant for length 
transition 


Number of coordinates 
X airfoil coordinates 


Y airfoil coordinates 
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of 


INPUT 
IPT (1) 
IPT (2) 
IRSTRT 
IGLMAX 
IPRINT 
RL 
XTRL 
XTRU 
ALPO 
GTR 
MPTS 
XP 


ye 


TABLE 6.2 


VALUES USED FOR THE CEBECI IBL CODE 


VALUE 
a 
1 


0 


101 
X airfoil coordinates 


Y airfoil coordinates 
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VII. PRESENTATION OF COMPUTATIONS 


Simulations of steady flow at Reynolds numbers ranging 
from 1 million to 15 million at zero degrees angle of 
attack, angle of attack studies for Reynolds numbers of 1.5 
Million and 6 million all over a NACA 0012 airfoil are 
presented in this chapter. 

Four aerodynamic factors are investigated: coefficient 
of lift, coefficient of pressure, velocity profiles at 
specified locations along the chord of the airfoil, and erin 
friction along the airfoil chord. In the steady flow cases 
results from the Sankar N-S code are compared with results 
from the Cebeci IBL code. Numerous studies have documented 
the validity of the Cebeci code in steady flow [Refs. 2,12, 
28,30}, however, the Cebeci IBL code only considers 
incompressible flow. Therefore, the Sankar N-S code was 
limited to the upper regions of what is normally considered 
incompressible flow, .3 Mach. Unless otherwise indicated, 
all cases are run at .3 Mach. 

The coefficient of lift as a function of angle of attack 
is plotted in Figure 7.1 for a NACA 0012 airfoil at a 
Reynolds number of 6 million. Figure 7.2 shows lift 


coefficients versus angle of attack for a Reynolds number of 


i> million. 
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Next, the correlation between the suction pressure 
coefficients for various Reynolds numbers and angles of 
attack were investigated. At zero degrees angle of attack 
the coefficient of pressure, Cp, of the upper surface is the 
same for the five Reynolds numbers presented: 1, 3, 6, 10, 
and 15 million, using both Cebeci's IBL code and Sankar's 
N-S code. Figure 7.3 plots this. Coefficients of pressure 
at angles of attack of 0O, 4, 8, and 12 degrees at a Reynolds 
number of 6 million and angles of attack of 0, 2, 4, and 6 
degrees for a Reynolds number of 1.5 million and a Mach of 
wi2 are shown in Figures 7.4 and 7.5 respectively. The 
conditions of Figure 7.5; 1.5 million Reynolds number, .12 
Mach, and 0, 2, 4, and 6 degrees angle of attack, were 
chosen to verify steady state conditions presented in Tang 
(Ref. 23]. The pressure coefficients were compared, 
however, no other results were DECeeneea by Tang. 

Skin friction and velocity profile information was then 
sought. The verification studies done by References 3, 7, 
14, and 15 did not address either. Therefore, the only 
comparisons made were with the Cebeci IBL code. Figure 7.6 
shows the coefficient of skin friction, Cg, as a function of 
alxrftolls chord. The Sankar N-S code is calculated as 
turbulent over the entire airfoil [Ref. 13]. The Cebeci IBL 


code allows specification of transition from laminar to 
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SKIN FRICTION 


_ Ts 





qe 


Ce vs x/c; Re = 6M, AOA = 0, N-S M 


Figure 7.6 


turbulent flow or allows the option to have the code compute 


transition [{Ref. 30]. 


Figure 7.6 plots Cre for the NACA 


0012 airfoil at a Reynolds number of 6 million at zero 


degrees angle of attack for eight cases, three using the 


Cebeci IBL code and five using the results of the Sankar N-S 


code: 


1) 
2) 


3) 


a) 
5) 
6) 
7) 
8) 


The 


Cebeci IBL Code 
turbulent flow; transition at .005c 
computed transition at .27c 
laminar flow; transition at .70c 

Sankar N-S code 
first grid point at .000l1c 
first grid point at .00005c 
first grid point at .00002c 
first grid point at .00001c 
first grid point at .000005c 
C-grid generates 41 =points ine the direction, 


clustering them near the wall and stretching them further 


from the wall. 


Specifying the first grid point's location 


from the wall determines where the grid clustering begins. 


This study was repeated for a Reynolds number of 3 


Million at zero degrees angle of attack in Figure 7.8. Six 


cases are presented, two using the Cebeci IBL code and four 


using the results of the Sankar N-S code: 
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Cebeci IBL code 
1) turbulent flow; transition at .005c 
2) computed transition at .44c 
Sankar N-S code 
3) first grid point at .000l1c 
4) first grid point at .00005c 
5) first grid point at .00001ic 
6) first grid point at .000005c 
Figure 7.8 presents a comparison of five cases for a 
Reynolds number of 15 million at zero degrees angle of 
attack: 
Cebeci IBL Code 
1) turbulent flow, transition at .005c 
2) specified transition at .70c 
Sankar N-S code 
3) first grid point at .0001c 
4) first grid point at .00005c 
5) first grid point at .00002c 
Velocity profiles are plotted in Figures 7.9 through 
7.17. Figures 7.9 through 7.12 show velocity profiles for 
the conditions of Figure 7.6: Reynolds number = 6 million, 
angle of attack = 0, N-S Mach = .3. Figure 7.9 is the IBL 
velocity profile computed when transition is specified at 
.005c. Figure 7.10 is the laminar velocity profile resulting 
from specifying the IBL transition at .70c. Figure 7.11 and 


7.12 are the velocity profiles computed when specifying the 
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N-S first grid point from the wall at .0001lc and .00002c 
respectively. 

Figures 7.13 through 7.15 are the velocity profiles 
associated with the conditions of Figure 7.8: Reynolds 
number = 15 million, angle of attack = 0, and N-S Mach = 
p30 Figure 7.13 is the IBL turbulent velocity profile, 
Figure 7.14 is the N-S turbulent velocity profile resulting 
from specifying the first grid point off of the wall as 
.0O001lc, and Figure 7.15 is the N-S velocity profile 
resulting from specifying the first grid point off of the 
wall as .000005c. 

Figures 7.16 and 7.17 are velocity profiles generated 
from a Reynolds number of 1 million at zero degrees angle of 
attack and .30 Mach number. Figure 7.16 is the IBL velocity 
profile resulting from specifying transition at .005c and 
Figure 7.17 is the N-S velocity profile resulting from 
specifying the first grid point at .000lc from the wall. 

The results of Figures 7.6 through 7.8 are replotted in 
Figures 7.18 through 7.20. However, only the turbulent 
cases and small grid mesh conditions are shown. Figure 7.18 
plots the IBL turbulent skin friction at zero degrees angle 
of attack for Reynolds numbers of 3, 6, and 15 million. 
Figure 7.19 plots the N-S skin friction at the same condi- 
tions as Figure 7.19 for the results of the case with the 
first grid point specified closest to the wall. Figure 7.20 


compares the IBL and N-S results from Figures 7.18 and 7.19. 
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Next, comparisons between the IBL and N-S skin friction 
for varying angles of attack were investigated for the 
conditions of Figure 7.4: Reynolds number = 6 million, Mach 
= .3, and Figure 7.5: Reynolds number = 1.5 million, Mach = 
.12. Figures 7.21 through 7.24 are for a Reynolds number of 
6 million and a Mach of .3. Figure 7.21 plots the IBL skin 
friction profile for angles of attack of 0, 4, 8, and 12 
degrees. Transition is specified at .005c. Figure 7.22 
plots the N-S skin friction for the same angles of attack 
and a first grid point off the of the wall of .000l1c. 
Figure 7.23 is also the N-S skin friction, but the first 
grid point off of the wall is specified at .00001lc. Figure 
7.24 compares the IBL and N-S skin friction profiles for 12 
degrees angle of attack. 

Figures 7.25 through 7.28 present skin friction profiles 
for a Reynolds number of 1.5 million and a Mach of .12. 
Figure 7.25 is the IBL turbulent skin friction (transition 
is specified at .005c). Figures 7.26 and 7.27 are the N-S 
skin friction profiles computed when the first grid point is 
specified at .000lc and .00005c, respectively. Figure 7.28 
compares the IBL and N-S skin friction profiles for 6 


degrees angle of attack. 
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VIII. RESULTS AND DISCUSSION 


The Sankar N-S code has been compared to the Cebeci IBL 
code. Integrated values that are not strongly influenced by 
viscosity, such as coefficients of lift and pressure, 
correlate well, as shown in Figures 7.1 through 7.5, and 
with the pressure coefficients presented in Tang [Ref. 23]. 
However, the viscosity influenced boundary layer values, 
such as the coefficients of skin friction and velocity 
profiles, are very sensitive to the type of flow and the 
grid size. This is evident in Figures 7.6 through 7.8. The 
IBL laminar skin frictions are much lower than the turbulent 
values. Also, the influence of the grid mesh on the ability 
of the N-S code to compute turbulent values is shown. When 
the first grid point off of the wall is chosen as .000l1c, 
the resulting skin frictions appear laminar for the 6 
million and 15 million cases (Figures 7.6 and 7.8). The 
lower the Reynolds number, the less sensitive the 
computations are to grid size, as can be seen by comparing 
the IBL and N-S skin friction values of a Reynolds number of 
3 million (Figure 7.7). In all cases, specifying the first 
grid point closer to the wall positions more points in the 
boundary layer. For Figure 7.6, when the first grid point 
from the wall is specified at .000lc, 9 points are located 


in the boundary layer at 10% chord, 16 points at 50% chord, 
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and 29 points at the trailing edge. In comparison, when the 
first grid point from the wall is specified at .000005c, 13 
points are located in the boundary layer at 10% chord, 19 
points at 50% chord, and 29 points at the trailing edge of 
the airfoil. 

When the influence on skin friction of subsequent grid 
refinements in the N-S code is negligible, the N-S skin 
frictions are substantially greater than the values computed 
by the IBL code when transition is specified at .005c. See 
Figure 7.20. It is of interest to note, however, the higher 
coefficients of skin friction evident when transition is 
delayed in the IBL code as in Figure 7.6. 

The velocity profiles, Figures 7.9 through 7.17, vary 
little in the turbulent cases, regardless of grid size or 
Reynolds number. The shape of the profiles varies between 
the IBL and N-S cases, but both exhibit turbulent velocity 
profiles. There is a discrepancy between the skin friction 
and velocity profile results computed by the N-S code when 
the first point off of the wall is .000lc. The coefficients 
of skin friction indicate laminar flow, whereas the velocity 
profiles are turbulent. A laminar velocity profile, Figure 
7.10, is calculated using the IBL code (Reynolds number of 6 
million, O degrees angle of attack). 

The influence of the grid mesh is again evident in the 
angle of attack studies, Figures 7.21 through 7.28. The 


coefficients of skin friction for the higher Reynolds number 
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of 6 million seem to not be dependent on angle of attack if 
the grid points are not located sufficiently close to the 
wall. This is not evident for the lower Reynolds number of 
1 million. The N-S skin friction values, when no longer 
influenced by grid size, are again higher than the fully 
turbulent values computed by the Cebeci IBL code, as seen in 
Figures 7.24 and 7.28. 

Varying the artificial viscosity in the Sankar N-S code 
did not vary the solution of the skin friction. However, 
all cases were run with the first grid point .0001ic from the 
wall, so no conclusion should be drawn. If the number of 
time steps was exceedingly large (8000 steps) for zero 
degrees angle of attack, the pressure and skin friction 
profiles began to indicate separated flow. 

Angle of attack studies using the Sankar N-S_ code 
required a large number of time steps (approximately 7000) 
before the coefficient of lift converged to realistic 
values, regardless of the first grid point off of the wall 
specified. The coefficient of friction values stabilized 
relatively quickly, within 2000 to 3000 time steps. Several 
computer runs at each angle of attack and Reynolds number 
were then required to determine when specification of the 
first grid point off of the wall no longer influenced the 
skin friction values computed. At this point, the results 


were assumed to be converged. 
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The IBL code took less than 16 iterations to converge, 
regardless of the Reynolds number chosen or the angle of 
attack. Very low Reynolds numbers would require more 


iterations, however. 
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IX. CONCLUSIONS AND RECOMMENDATIONS 


Both the Sankar Navier-Stokes Code and the Cebeci 
Interactive Boundary Layer code show reasonable results for 
high Reynolds number, incompressible steady flows over a 
NACA 0012 airfoil. The discrepancy in skin frictions should 
be resolved and the influence of dissipation on the skin 
friction investigated. These steady-state results should 
then be extended to other airfoils and low Reynolds flows to 
determine the effect of increased viscosity on the codes. 
The effect of transition on the velocity profiles and the 
skin friction is modeled in the Cebeci IBL code; however, 
the lack of a transition mode: and a smearing function in 
the Sankar N-S code limits its ability to model most 
experimental flows, especially at low Reynolds numbers. 
Also of importance is the wake influence, which has not been 
addressed in this report, and the growth of the boundary 
layer along the airfoil. Since the velocity profile results 
are presented non-dimensionally, this trend is not evident. 
The difference in profile shapes generated from the N-S code 
and the IBL code could be resolved better if the actual 
boundary layer profile, rather than a non-dimensionalized 
profile, is used. 

Neither code was extensively compared to experimental 


results for skin frictions or velocity profiles. This is, 
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of course, a criterion in determining the accuracy of the 
codes. The cost and time considerations associated with 
running the codes indicate that, at the present, for 
steady-state runs, the IBL code is more effective. The N-S 
runs were submitted to the NASA X-MP Cray, where waiting 
times of 24 hours before the program was executed was 
common. Execution times ran between 10 and 40 minutes, at a 
cost of $1000 per hour. In comparison, the IBL program 
normally was completed in 10 to 30 minutes on the Naval 
Postgraduate School IBM mainframe, for less than $50. 

The N-S code is a very effective tool for calculating 
dynamic stall characteristics. Its capability as a general 
flow solver, however, is limited in comparison to the IBL 


technique in terms of cost and time constraints. 


102 


10. 


LIST OF REFERENCES 


McCroskey, W.J. and McAlister, K.W., et al., "Dynamic 
Stall on Advanced Airfoil Sections," J. American 
Helicopter Society, V. 26, No. 3, pp. 40-50, July 1981. 


Cebeci, T., "Airfoils with Separation and the Resulting 
Wakes," J. Fluid Mechanics, V. 163, pp. 323-347, 1986. 


American Institute of Aeronautics and Astronautics Paper 
85-0129, "Numerical Solution of Unsteady Viscous Flow 
Past Rotor Sections," by N.L. Sankar and W. Tang, 
January 1985. 


Beam, R.M., and Warming, R.F., "An Implicit Factored 
Scheme for the Compressible Naviar-Stokes Equations," 
AIAA Journal, V. 16, No. 4, pp. 393-402, April 1978. 


Steger, J.L., "Implicit Finite-Difference Simulation of 
Flow about Arbitrary Two-Dimensional Geometries," ATIAA 
Journal, V. 16, No. 7, pp. 679-686, July 1978. 


American Institute of Aeronautics and Astronautics Paper 
87-0495, “Boundary Layer Measurements on an Airfoil at 
Low Reynolds Numbers," by M. Brandal and T.J. Mueller, 
January 1987. 


Marchman, J.F., "Aerodynamic Testing at Low Reynolds 
Numbers," J. Aircraft, V. 24, No. 2, pp. 107-114, 
February 1987. 


National Aeronautics and Space Administration Technical 
Memorandum 84245, "An Experimental Study of Dynamic 
Stall on Advanced Airfoil Sections," by W.J. McCroskey, 
K.W. McAlister, L.W. Carr, and S.L. Pucce, July 1982. 


American Institute of Aeronautics and Astronautics Paper 
87-0242, "A Study of the Laminar Separation Bubble on an 
Airfoil at Low Reynolds Numbers Using Flow Visualization 
Techniques," by G.S. Schmidt and T.J. Mueller, January 
1987. 


McCroskey, W.J., Carr, L.W., and McAlister, K.W., 


"Dynamic Stall Experiments on Oscillating Airfoils," 
AIAA Journal, V. 14, No. 1, pp. 57-63, January 1976. 


103 


ieee 


Ze 


13. 


14. 


15. 


16. 


aly. 


18. 


19. 


20. 


PAN 


PAG) & 


23 ¢ 


Cebeci, T., Editor, Numerical and Physical Aspects of 
Aerodynamic Flows, Vols. I-III, Springer-Verlag, New 


York, 1982, 1984, 1985. 


Subroto, P.H., Viscous/Inviscid Interaction Analyses of 


the Aerodynamic Performance of the NACA 65-213 Airfoil, 
Master's Thesis, Naval Postgraduate School, Monterey, 


California, March 1987. 


Valdes, James, Dynamic Stall Calculations Using a 
Navier-Stokes Solver, Master's Thesis, Naval 


Postgraduate School, Monterey, California, December 
1986. 


Daily, J.W. and Harleman, D.R., Fluid Dynamics, Addison- 
Wesley Publishing Company, Inc., 1966. 


Shapiro, A.H., The Dynaics and Thermodynamics of 
Compressible Fluid Flow, Vol. 1, John Wiley and Sons, 


DNC. melo SD 3 


Schlichting, H., Boundary Layer Theory, 7th ed., McGraw- 
Hill, Book Co., 1979. 


Moran, es? An Introduction to Theoretical and 


Computational Aerodynamics, John Wiley and Sons, Inc., 
1984. 


Patel, V.C., Rodi, W., and Scheuerer, G., "Turbulence 
Models for Near-Wall and Low Reynolds Number Flows: A 
Review," AIAA Journal, V. 23, No. 9, pp. 1308-1319. 


Cebeci, T., Chang, K.C., Li, C., and Whitelaw, J.H., 
"Turbulence Models for Wall Boundary Layers," AIAA 
Journal, V. 23, No. 3, pp. 359-361, March 1986. 


Cebeci, T. and Smith, A., Analysis of Turbulent Boundary 
Layers, Academic Press, 1974. 


York, B. and Knight, D., "Calculation of Two-Dimensional 
Turbulent Boundary Layers Using the Baldwin-Lomax 
Model," AIAA Journal, V. 23, No. 12, pp. 1849-1851, 
December 1985. 


American Institute of Aeronautics and Astronautics Paper 
81-1289, "Dynamic Stall of NACA 0011 Airfoil in 
Turbulent Flow--Numerical Study," by Y. Tasner and N.L. 
Sankar, June 1981. 


Tang, W., Numerical Solutions of Unsteady Flow Past 
Rotor Sections, Ph.D. Thesis, Georgia Institute of 


Technology, Atlanta, Georgia, August 1986. 


104 


24. 


2's 


26. 


Oa ae 


28. 


29. 


30. 


Bl. 


Jameson, A., "Iterative Solution of Transonic Flows over 
Airfoils and Wings, Including Flows at Mach 1," 


Communications on Pure and Applied Mathematics, V. 
XXVil, pp. 283-309, 1974. 


Howard, R., Class notes from High Speed Aerodynamics, 
Naval Postgraduate School, Monterey, California, Winter 
1987. 


Murman, E. and Cole, J., "Calculation of Plane Steady 
Transonic Flow," AIAA Journal, V. 9, No. 1, pp. 114-121, 
January 1971. 


Cebeci, T., Wang, G.S., Chang, K.C., and Choi, Jd., 
"Recent Developments in the Calculation of Flow over Low 
Reynolds Number Airfoils," paper presented at _ the 
International Conference on Aerodynamics at Low Reynolds 
Numbers, held at the Royal Aeronautical Society, London, 
England, 16-17 October 1986. 


Keller, H.B., "Accurate Difference Methods for Two-Point 
Boundary Value Problems," SIAM J. Numerical Analysis, V. 
11, No. 2, pp. 305-320, April 1974. 


Paik, S.W., Application of the Viscous  Inviscid 
Interaction Methods to the Analysis of Airfoil Flows at 


Low Reynolds Numbers, Master's Thesis, Naval 
Postgraduate School, Monterey, California, December 


1986. 


Anderson, D.A., Tannehill, J.C. and Pletcher, R.H., 


Computational Fluid Mechanics and Heat Transfer, McGraw- 
Hill (Hemisphere Publishing Co.), 1984. 


Halsey, N.D., "Potential Flow Analysis of Multielement 


Airfoils Using Conformal Mapping," AIAA Journal, V. 17, 
No. 12, pp. 1281-1288, December 1979. 


105 


INITIAL DISTRIBUTION LIST 
No. Copies 


Defense Technical Information Center 2 
Cameron Station 
Alexandria, Virginia 22304-6145 


Library, Code 0142 2 
Naval Postgraduate School 
Monterey, California 93943-5002 


Chairman, Code 67 5 
Department of Aeronautical Engineering 

Naval Postgraduate School 

Monterey, California 93943-5004 


Lawrence W. Carr AL 
Fluid Mechanics Laboratory 

Mail Stop 260-1 

NASA Ames Research Center 

Moffett Field, California 94035 


N.L. Sankar 1 
School of Aerospace Engineering 

Georgia Institute of Technology 

Atlanta, Georgia 30332 


Tuncer Cebeci al 
Director, Center for Aerodynamic Research 
California State University 

Long Beach, California 90840 


Campbell Henderson 2 
Head, Aircraft Performance Analysis 
Branch 
Code 6051 
Naval Air Development Center 
Warminster, Pennsylvania 18974-5000 


Lisa J. Cowles 10 
Code 6051 | 

Naval Air Development Center 

Warminster, Pennsylvania 18974-5000 


106 

















thesC756965 
High Reynolds number, low Mach number, 


ON) 


3 2768 000 77208 1 


DUDLEY KNOX LIBRARY 


















































































































































Siisg ' y ° 
é . ' 
1 
. . ' 
1 "4 ? ' 
‘ 
. e on 
. 4 
. t ‘ 
° ; 
' 
' . 
. “1 , an ' 
. ' ’ 
8 . 
s te ma . 
° ' . 
. ‘ ’ ' 
é ' . 
: } 4 
2 C . . ie 
a al ® ' 
* © Beet Gpahte 2 BK Rts Pe hee 7 2 ‘ ' . ’ or 8 
FP IOS AOS 6d ph. hal ibis Boke , tt y Phe . ' 
C2 Ce 4 oe pee Pith A ae hd ‘8 $45 ‘ = ' '‘ LJ 
ats PM Be, 64 6 ayes month ods om J ' . 
“ot homme: das aD 0 ate Ae mile) tous ’ ¢ . ‘ . ' 
Pht Bk anes SPR AD Dh 50 eo Op e ' ‘ ’ 
2.2% (nar Fgh 4 ? hr Dee Ae “eo © g's ’ ’ ‘ 
P28 ee tetematne Oh. Orhee .- . ' 
#F.. Ch2atirn me "Praise de Sai a x ; 
2 ne atau Loe ee "* Bo soy ' 
“* Sel @ Beto tnn oy, watt of s 
- Dw os Pom. aw s 
22 PS ar ' f a 
© 7 4 e ie s ' 
Sut, ‘ . 
ah A 10 pony: ¢ ' ® 
4 @ a eo , 
“at op ah ' . ‘ 
- 7 ‘ ' ' . ’ 7 
BOE te. ane ge’, , - - ' , e 
oe U t ’ . ' 
oa) ' ; ' ‘ 
0 all oF 4a "’ ’ 
Pe as, ' ’ s m4 ’ 
Pink Catt » : : 
at A ar Rte ae oe E u 
i nietentitme atin te .! 
ODO ee al nee : : 
tarot i, 7 7 
= Got +4 Of tn Mae 5 . 
ome f¢ iy ’ ' r ' 
FP FS lt wth to « ha ’ ‘ F 5 
MOOD ttt ho ca , ® 
a ate le a oy ’ ¢ ' 
POM, atoll ay amy op oo L 2 , 
‘ ’ 
‘ ’ ' 
sa ‘ ' 4 
L 5 ft ri¢ . ' 
AA ep ' ’ ' 1 . ‘ 1 1 * 
ar) . o. 
er | . ’ oo 
’ ‘ ' 
5 ' : ° ' 
' ‘ ' to. 
’ a ’ ' ' = ' 
oe ' ' ® ‘ ‘ 
“3 ' ‘ ' 
‘ ‘ . 
' ' 
* te 4 . 
‘ ou ’ 
e ry 
’ 
: ' 1 
’ , , 4 ' ‘ 1 ’ 
' ‘ ' ' 
s of ' 
; 7 ' 
‘ ' ' ' 
' ’ . e wos 
: ’ ‘ ' 
; cy ‘ 
’ 4 , U ’ 
‘ ' . ou ' . 
a ‘ . 7 ' 
' ' * ' 
' ‘ ' 
' ' 
' ' ' 
. * e . 
ee . an | 5 cd ‘ id ' 
r ’ ’ ‘ 
= ‘ ‘ ° . ' 
. . F if ' 
, ' - 
' ' ' J] 1 i 
’ = » . 
é ' ‘ J ' . ’ 
es B ’ 
: . : ' ‘ 
a? . 8 ' 
s ’ > ° ~ ' 7 ' . ‘ ' 
' ' . ’ 7 ' 
‘ , ) ; ’ 
5 ba ‘ ‘ ' ‘ . 
' 
A : : : 
' ' 
; ts - . : 
5 s) ’ , 
‘ : é ‘ 
: , . 8 ‘ 
® 4 : ' 
: 
, . ' "4 . ' 1? 
My 
' P 7 4 ' 
' 1 ' ' 
’ 
: m8 ' . 
1 a ’ ' 
r os ; to. ' 
Ps ’ ' ' ' ' ' 
~ : U é 
' ’ : é 
' ’ ° 
ae ' f 6 
: ’ 
’ ? ‘ ‘ 
’ e 
e ? e 
e s . 
e ' ' 
’ 2 " ’ 
« . 
' ' ’ ' ° . # 
; : 
Le 4 . bd os 
’ 4 ' . 
: to? 8 roy ' ' " ‘ : 
' & . s ‘ s ‘ st 
*s OF Jr Bs , vous : ' 
' ' ’ 4 Co Se) : . 
. ' . ' s ' 1 
. e" . ‘ * . ' 
ase ¢e,0' 4 we ’ % Fe ‘ , 
a Rrrerig y ar) wise e “1g ' ' ' 
© se : e » ' +8 e = 3 a c ry Pi Lat 7 1 e ; 
(em « “S* $ « Luss ' ee é ¢ ‘ 
4 9° ewe, ye he Ie FO wte.* we bet ny ee a a, a = . ’ a ’ ' . e ' 
P bth bitecte Yorke te TT OX oman wy er ‘ » %~ ge 4 oe 6te . ee ’ 7) o ' ' 
a) oO S Sek, ree ew 4 eee 89s OF OEE) ore ae te Risers ‘ for L ’ t ' ' , ‘ 
“=e > 26 ne! oe © Smeg mee Oa ce 4 ts =) 6 e% ~ '>"R se Ld ' es . s ‘ ’ ' . . 
|e ot, SC em a Le MH Peredsa, Oe! caee és pute e 1 e es e 4 a ’ 4 1 os 
FOOT ee ee wo ay See eed, Laie 1 Ys . i | ‘e ss 7 uh ae | ' ' ‘ ; ' a ' ' 
Fore we Ps me @ eet pene ‘ > ' be eo . *Be . ee U a¢ r ' . 
ed. ee ’ ‘4 ’ + t . em ' Se , b é - ‘ ; 
28 SS we ew ve eS wate 19, > en ’ t er a ’ ° ss ‘ a ' ‘ 
mE! wy Ouse, wags me , fe Dever ' 5°@ omy . ry ‘ 
Te Me one wm ‘ tT | in a 4 sere - . er ae | ' e ' aos ‘ ; 
ONE OM te Re ey: a See e tr @ ee be rt ’ pat ‘ . . e . b] 
. 6 SETS She OE DHE: woe hE Op recess , es eirde @ie *04 © b w% ‘ : ’ ' ’ * ' ' ' ' ' ' 
oye Pet yee ry we - & woe >) 6m e 4 e! 7? Q ’ * ' Q 
ee ele 4 1 ¢@ ' . ' ' » 
] ‘es i a Sa Sree e s ‘ 
4 West fare « * win Da ar ’ ’ ° e te 
“ y os . ' Ls ya ‘ aos , ' ' e ' 
Stee eis i a vie o 4 at e , “! 
s ' ~o ' Phe gra ' ‘ ' ’ ‘ ’ . 
»é a * | | e'e@ ef . e 7 4 « ' ' ' 
**. e606 ™=% Ge, Btls OPoass oe i oY SVerea for e @e& oF a ° , 4 * er é . 
oy 5 Oa Pee % Pes FOr te ape s eo at a as ee ‘ ra 4 5 
’ 7 % 56 Siggy oT aS See 6 ay are ‘ ry ° . ore ° ° ‘ 
ae for ay) Ose Oy ag a %& fe ae 4 @& tere ere I an oe ' : , ® ‘ 
| TS Ue Oe haw wey Dace ta ce © , ee 6% af is Ge ’ 7 
PPR rE ea thee weg ’ ea % 8 be a ‘ >i Bese Sas ' ’ ‘ ® ' 
$60 6ace Sr ace t@ oe yas re ase 7 oy r ’ r tee ’ 
POSE Pure hy Stes (Pow, G\oe & weep « ‘ ry .Y ' es ge ‘ 
 0°R © vera On °@ Pk case ng Swe -~ see a , : . 
S29 e eer “S'S Pal «ee ‘ “ 4 ® ‘ ’ e a. ’ ' 
88's ares ees Mee e && Ti] r) p ' . 24 
COPS PEO sere. y ‘ if ‘ » ‘ ’ ' ' ’ ' ? 
. 1 ' 
) ' t ' 
@e '@-e 6 ww ee ' ; ’ . 
'*@ bea eo Be wer > . r& \: s , ‘ 1 
8: em sat see ‘ ee ' ¢ ; 7 a ' on ' ope ’ 
7,0 Keane , s a, ‘4 ' inn 
' ‘ ; ‘ ’ ’ 
' ’ ' ' s ‘ ‘ ’ 
‘ é ‘ Uy é ‘ 
wwe Cee wert ' e ees Fy a ’ ' ® ' 
Severe taew\"¢ » i & ’ ' ‘ 
 @e Tee boom O@t Rrra ' ge ’ ' ‘ e e 4 ’ 
STOP: wPa ce o& "BL PPevs u (hes ' ' ; ' ‘ 
pe FO se OP ie * bean ’ Py ‘ ‘ 4 
6 Se aelyy 90 8s Bh oPetecs hee tL 4 ee s ‘ 
2 70° 9 weg a: » 'e q ' ‘ 
CFO"R CFS Ove ws, ‘ 4 , ’ 
OD FO7R 101 Os Oe ; ’ 4 4 . Z ‘ 
‘Sphe OH we we Fate MHEG Mew hes's "ew eases as t ’ ‘ e ’ ’ s ° ' 
ve “h WeeTet., @ ' t s L 
‘ a , 
q . | es ‘ ¢ ® 
‘ : ’ ‘ 
? ’ ‘ 
es es. et , e 
ry & Fy ' ore ' 
‘uy ry ‘ ' 
' ' 8 * 6 
wa Os "6a © % ay ¥ ian LJ ‘ ' 5 
b thet dd le dt TTT » a r ‘ , ‘ 
tr at 9 5 .? ‘ ' 
7» S't @ ‘ a P TT a? ' ' 
r 7°S™ we Foe Cops T'Oe og e ‘ os t ‘ ' 
> Ole were ees ? ry ‘ ‘ 4 ‘ s ‘ rs 
were :ee c B ‘ 
, oPhee i t ' ' ¢ ee , t ; 
#4 “ee ' 
MEO Oru ; ' 
a ‘ ‘ 
wre rres 8 8 . + 
4 ®vre ' ® ‘ 
roe 7 ’ 
eed Att 
oe fy 
Carpep 
om ‘ 
CFP eee wwe i) 
GOR —oOV ID we, CPE Ms Bg bite Sk oe Te Er ree 
SEOs O tayre edie, by ree om ud T7Weaos as. 
steht ty he 2 to te EL ee 





